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ABSTRACT 

Aims. Three dimensional explicit hydrodynamic codes based on spherical polar coordinates using a single spherical polar grid suffer 
from a severe restriction of the time step size due to the convergence of grid lines near the poles of the coordinate system. More 
importantly, numerical artifacts are encountered at the symmetry axis of the grid where boundary conditions have to be imposed that 
flaw the flow near the axis. The first problem can be eased and the second one avoided by applying an overlapping grid technique. 
Methods. A type of overlapping grid in spherical coordinates is adopted. This so called "Yin- Yang" grid is a two-patch overset 
grid proposed by Kageyama and Sato for geophysical simulations. Its two grid patches contain only the low-latitude regions of the 
usual spherical polar grid and are combined together in a simple manner. This property of the Yin- Yang grid greatly simplifies its 
implementation into a 3D code already employing spherical polar coordinates. It further allows for a much larger time step in 3D 
simulations using high angular resolution (:£ 1°) than that required in 3D simulations using a regular spherical grid with the same 
angular resolution. 

Results. The Yin- Yang grid is successfully implemented into a 3D version of the explicit Eulerian grid-based code PROMETHEUS 
including self-gravity. The modified code successfully passed several standard hydrodynamic tests producing results which are in 
very good agreement with analytic solutions. Moreover, the solutions obtained with the Yin- Yang grid exhibit no peculiar behaviour 
at the boundary between the two grid patches. The code has also been successfully used to model astrophysically relevant situations, 
namely equilibrium polytropes, a Taylor-Sedov explosion, and Rayleigh-Taylor instabilities. According to our results, the usage of 
the Yin- Yang grid greatly enhances the suitability and efficiency of 3D explicit Eulerian codes based on spherical polar coordinates 
for astrophysical flows. 

Key words. Methods: numerical - Hydrodynamics - Gravitation - Supernovae: general 



1. Introduction 

Three dimensional hydrodynamic simulations employing a sin- 
gle spherical polar grid are computationally expensive be- 
cause of the convergence of grid lines towards the north and 
south pole. The converging grid lines imply a severe restric- 
tion of the time step size for any hydrodynamic code using 
explicit time discretization due to the CFL condition. This 
so-called "pole problem" bothers astrophysicists when simu- 
lating self-gravitating flow in three dimensions (e.g., convec- 
tion in stars, or stellar explosions) where the spherical coor- 
dinate system is often preferable. In particular, simulations of 
core-collapse supernovae are a problem with which astrophysi- 
cists have been struggling. While observations show clear ev- 
idence of asymmetric (3D) complex structures in supernova 
ejecta, numerical simulations, in most cases, are carried out 
only in two spatial dimensions assuming axisymmetry (e.g., 
Blondin & Mezzacappall2006t IScheck et al Jl2006t: lOhnishi etall 
2007). Three dimensional core-collapse supernova simulations 
are rare (e.g., Janka et al. 2005; Mezzacappa et al. 2006; Scheck 
I2006t [iwakami et al.ll2008l)~ In addition to the severe restriction 
of the time step size, boundary conditions that have to be im- 
posed at the symmetry axis 8 e [0, n] flaw the simulations near 
the axis by causing undesired numerical artifacts in 2D axisym- 
metri c simulations, as e.g., jet-like flow features ( Kifo nidis et al.l 
120031 In 3D simulations, the axis represents a coordinate sin- 
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gularity that almost unavoidably will leave its mark on the flow 
near or across the axis. 

There have been attempts to construct a new type of grid 
which is able to ease the pole problem. However, it is not pos- 
sible to construct a single grid patch that can cover the en- 
tire surface of a sphere, is orthogonal, and at the same time 
does not contain any coordinate singularity except at the ori- 
gin. Therefore, multi-patch grid and overlapping (or overset) 
grid approaches are employed. They are widely used in the field 
of computational fluid dynamics where complex grid structures 
are common. For flows possessing an appro ximate global spher- 
ical symmetry, the "cubed sphere" grid dRonchi et alJ [l996) 
has been develope d and is currently applied to several astro- 
physical problems dKoldoba et al. 2002t iRomanova et al . 2003; 
IZink et alJl2008t iFragile et al.ll2009L e.g.). It is an overset grid 
consisting of six identical patches covering a solid angle of 4n 
steradians. The "Yin- Yang" grid has the latter property, too, but 
up to now it has not been used in astrophysical applications. 

Th e Yin- Yang grid was introduced by iKagevama & Satol 
(2004). It consists of two overlapping grid patches named "Yin" 
and "Yang" grid. In comparison with other types of overset grids 
in spherical geometry, the Yin- Yang grid geometry is simple, 
as both the Yin and the Yang grid consist of a part of a usual 
spherical polar grid. The transformation of coordinates and vec- 
tor components between the two patches is straightforward and 
symmetric, thus allowing for an easy and straightforward imple- 
mentation of the grid into a 3D code already employing spher- 
ical polar coordinates. The Yin- Yang grid is successfully used 
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on massively parallel supercomputers in the field of geophysi- 
cal science for simulations of mantle convection and the geo- 
dynamo. In these applications the thermal convection equation 
and the magnetohydrodynamic (MHD) equations are solved on 
the Yin- Yang grid using a second-order accurate finite difference 
method. Here, we also adopt the Yin- Yang grid, and use it for as- 
trophysically relevant (finite-volume) hydrodynamic simulations 
for the first time. 

The paper is structured as follows. In section 2, we describe 
the basics of the Yin- Yang grid configuration including the trans- 
formations of coordinates and vectors between the Yin and Yang 
grid patches. In section 3, we provide the details of the imple- 
mentation of the Yin- Yang grid into the PROMETHEUS hydro- 
dynamic code, and also discuss the resulting necessary modifica- 
tions of its 3D gravity solver that is based on spherical harmon- 
ics. In section 4, we present the results of the test calculations 
we have performed including a test with self-gravity. In section 
5, we discuss the conservation problem arising when applying 
the Yin- Yang grid. Then we report on the efficiency and per- 
formance gain obtained with the Yin- Yang grid compared to a 
spherical polar grid in section 6. Finally, we give the conclusions 
from our study in section 7. 



2. Yin-Yang Grid 

The Yin- Yang grid configuration is shown in Fig.Q] Both the Yin 
and the Yang grid are simply a part of a usual spherical polar grid 
and are identical in geometry. The angular domain of each grid 
patch is given by 
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where 8 and tp are the colatitude and azimuth, respectively. Note 
that it is necessary to add at least one extra buffer grid zone to 
both angular directions in order to ensure an appropriate overlap 
of the grids. The angular width 6 of this buffer zone depends on 
the grid resolution, i.e., 6 = AO = A<f>, where for simplicity we 
assumed equal angular spacing in 9- and 0-direction. The angu- 
lar domain is hereby extended by 25 in both angular directions. 
The Yin and Yang grid are patched together in a specific man- 
ner forming a spherical shell with a small overlapping region 
covering approximately 6% of a sphere's surface. Stacking up 
Yin- Yang shells in radial direction results in a 3D grid that is 
identical to the usual spherical polar grid in radial direction. It 
is obvious that, unlike in the case of the spherical polar grid, the 
problematic high latitude sections of the sphere are avoided, and 
the angular zoning is almost equidistant. 
The Cartesian coordinates 
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corresponding to the Yin grid, denoted by a superscript («), and 
the Cartesian coordinates 

(x (e \y {e \z (e) ) = f/sin0 w cos0 w , rsin(9 (e) sin0 (e) , rcos0 w ) 

(3) 

corresponding to the Yang grid, denoted by a superscript (e), are 
related to each other through the transformation 
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Fig. 1. An axis-free Yin- Yang grid configuration plotted on a 
spherical surface. Both the Yin (red) and Yang (blue) grid are 
the low latitude part of the normal spherical polar grid and are 
identical in geometry. The Yang grid is obtained from the Yin 
grid by two rotations, and vice versa. 



where 
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This Yin- Yang coordinate transformation can also be considered 
as two subsequent rotations. Accordingly, the transformation 
matrix M can be written as /? A (90 o )/? z (180°), where R x and R z 
are the transformation matrices of rotations by 90° around the x- 
axis and by 1 80° around the z-axis in counterclockwise direction, 
respectively. For the inverse transformation matrix M~ l = M 
holds. 

The relation between the spherical coordinates of the Yin and 
Yang grid patches can be derived directly from the transforma- 
tion matrix M. Because the Yin- Yang coordinate transformation 
involves only rotations, it implies that the radial coordinate is 
identical on the Yin and the Yang grid. The angular coordinates 
transform as 

6> (£,) = arccos(sin6» (,,) sin^ ( " ) ), 

cos ef n) 
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Note that the inverse transformation has the same form as © 
and (0 but exchanging the (grid) superscripts. 

Vector components in spherical coordinates transform ac- 
cording to 
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is the vector transformation matrix. When switching (grid) su- 
perscripts (e) and («) in matrix P, the inverse vector transfor- 
mation matrix is obtained. For a detailed deri vation of the trans- 
forma tion matrix P, we refer to section 3 of Kagevama & Satol 
(2004). Note that the vector transformation matrix P is singular 
at sin# w = 0, but this singular point is rectifiable. In practice, 
one can always decompose vectors into their Cartesian compo- 
nents and perform the corresponding transformation. 



3. Implementation 

We have implemented the Yin- Yang grid into our explicit finite- 
volume Eulerian hydrodynamics code, PROMETHEUS, which 
integrates the equations of multidimensiona l hydrodynamics us- 
ing th e piecewise parabolic method (PPM: ICollela & Woodward! 
1984) and dimensional splitting. The code also includes a 
Poisson solver based on spherical harmonics to handle self- 
gravity. 



3.1. Hydrodynamics solver 

Firstly, the Yin- Yang grid needs to be constructed. Since both 
the Yin and the Yang grid are part of a spherical polar grid an 
analogous spatial discretization in angular direction can be used. 
For example, the 9 and <p coordinates of the zone center of an 
angular zone (J, k) of a Yin- Yang grid, having Ng zones in 6- 
direction and zones in 0-direction, are given by 



<t>k = (, 
where 
A9 = 

A(f> = 



A9 
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for 1 < j < No, 
for 1 < k < Nj, 
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are the respective angular grid spacings. 

The range of values for the colatitude 6 and the azimuth angle 
d> are as given in ([1]), and for simplicity we set A6 = Ad>. In radial 
direction no modification is required. The geometric property of 
the Yin- Yang grid allows us to make use of the coordinate arrays 
r,, 0j, and fa twice by enforcing the same grid resolution for 
both grid patches. This approach avoids doubling the coordinate 
arrays. 

Only simple modifications are needed concerning the data 
and program structure. Arrays with three spatial indices, e.g., i, 
j, and k, need an extra grid index, say, I. For example, the ar- 
ray for the density field will be p(i, j, k, I) instead of p(i, j, k). As 
a consequence any triple loop running over indices i, j, and k 
in the program becomes a fourfold loop over i, j, k, and / in- 
stead. Otherwise, the Yin- Yang grid allows one to exploit with- 
out any further modification any already implemented finite- 
volume scheme in spherical coordinates to solve the equations 
of hydrodynamics. 

Different from the spherical polar grid, the Yin- Yang grid 
requires no boundary conditions in angular directions. Each 
grid patch communicates with its neighboring patch using in- 
formation from ghost zones that is obtained by interpolation of 
data between internal grid zones of the neighboring grid patch. 
Interpolation is only required in the two angular coordinates 




Fig. 2. A Mercator projection of an overlap region of the Yin- 
Yang grid. In case of bi-linear interpolation, four neighboring 
values of the underlying grid (red) will be used to determine the 
zone-centered value of a ghost zone in the grid on top (blue). 
The interpolation coefficients are determined by the relative dis- 
tances, denoted by black lines, between the interpolation point 
(diamond) and the four neighboring points (crosses). 



as the radial part of the Yin- Yang grid is identical to that of a 
spherical polar grid. It is straightforward to determine the corre- 
sponding interpolation coefficients. The mapping of vector quan- 
tities between the Yin and Yang grid patches requires an addi- 
tional step. After interpolating the vector components they must 
be transformed according to the transformation given in Eq. ^ 
from the Yin to the Yang angular coordinate system, and vice 
versa. 

We tested two interpolation procedures. In the first one all 
primitive state variables (density, velocity, energy, pressure, tem- 
perature, abundances) are interpolated ignoring the resulting 
small thermodynamic inconsistencies. In the second procedure, 
we only interpolate the conserved quantities (density, momen- 
tum, total energy, and abundances), and compute the velocity 
and the remaining thermodynamic state variables consistently 
via the equation of state. Both procedures produce very similar 
results which differ at the level of the discretization errors. As 
the second procedure is more consistent we use it as the stan- 
dard one in our code. 

An example of overlapping situations which are encountered 
when using a Yin- Yang grid is shown in Fig. [2] For simplicity, 
we use bi-linear interpolation in order to prevent unwanted os- 
cillation. Because the grid patches are fixed in both angular di- 
rections the interpolation coefficients for each ghost zone need to 
be calculated only once per simulation at the initialization step. 
After initialization, the coefficient map is stored in an array for 
later usage. Moreover, the symmetry property of the Yin- Yang 
transformation allows one to make use of the interpolation coef- 
ficients twice for both grids. 

Because the Yin- Yang grid is an overlapping grid integral 
quantities such as the total mass or total energy on the compu- 
tational domain cannot be obtained by just summing local quan- 
tities from every grid cells. Doing so will result in counting the 
contributions in the overlapping region twice. To circumvent this 
problem, weights are given to each grid zone during the summa- 
tion. Suppose a grid zone has an overlapping volume fraction a 
the cell will receive a weight w = 1.0 - 0.5a. Zones in the non- 
overlapping region receive the weight of 1 .0, i.e., the entire zone 
contributes to the integral while, on the other hand, zones that are 
fully contained within the overlapping region have a weight 0.5. 
The volume fraction a does not depend on the radial coordinate 
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and can be thought of as an area fraction since the grid patches 
are not offset in radial direction. Prior to the area integration, 
one needs to determine for each zone interface of the underlying 
grid the points where the interface is intersected by the boundary 
lines of the other grid, e.g., points on the Yin grid intersected by 
the boundary lines of the Yang grid. The intersection points can 
be determined using the Yin- Yang coordinate transformation in 
(0 and O, respectively. The integration in the overlapping area 
is then carried out using the trapezoid al method. This procedure 
is also described in Peng et al] (120061) . Once the area or volume 
fraction a is calculated, the weights for each cell are obtained 
easily. Note that these weights need to be calculated only at the 
initialization step, and are stored for later usage in a coefficient 
map w(j, k), where j and k are the indices referring to the 6 and 
coordinates, respectively. The coefficient map can be applied to 
both grids without any modification. Using the above described 
approach, the volume or surface area of the grids can be calcu- 
lated with an accuracy up to machine precision. 

3.2. Gravity solver 

The 3D Newtonian gravitational potential is computed from 
Poisson's equation in its integral f orm using an expansion int o 
spherical harmonics as described in Mi iller & Steinmetzld 19951) . 
Because the algorithm of these authors is based on a (single) 
spherical polar grid the density on the Yin- Yang sphere has to 
be interpolated onto an auxiliary spherical polar grid. The inter- 
polation used is first-order accurate, and due to the simplicity 
of the Yin- Yang grid configuration has to be performed only in 
the two angular dimensions. Concerning the resolution of the 
auxiliary grid, it is natural to employ the same grid resolution 
as that used for the Yin- Yang grid in all three spatial dimen- 
sions. The orientation of the auxiliary grid can be chosen freely 
in principle. However, it is convenient to align it with one of 
the two grid patches (the Yin-grid in our case). Once the den- 
sity is interpolated onto the auxiliary spheri cal grid we compute 
the gr avitational potential, as suggested by Mull er & Steinmeta 
( 1995|), at zone interfaces instead of at zone centers on both the 
Yin and Yang grid. The gravitational acceleration at zone cen- 
ters can then be obtained by central differencing the potential. 
Note that the interpolation coefficients for the density need to be 
calculated only once per simulation, because both the auxiliary 
grid and the Yin- Yang grid are fixed in angular directions. In ad- 
dition, all angular weights, Legendre polynomials, and their in- 
tegrals required for the calculation of the gravitational potential 
are stored after the initialization step for later usage. 

It is also possible to directly calculate the gravitational ac- 
celeration at zone centers. The gravitational potential i s given by 
(see Eqs. (5), (6) and (7) in iMuller & Steinmetzl d 19951) ) 

°° A 1 /I 

<D(r, 9, 0) = -G J] [ Z Y "" (e > ^ h7T C '» + r ' D ""^ 

1=0 m=—l * 
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where Y 1 '" and Y ,m * are the spherical harmonics and their com- 
plex conjugates, p is the density, and d Q. = sin 8 dd dip. The grav- 



itational acceleration in radial direction is then 
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Writing the radial derivative in Eq. (fTTI i as 
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and noticing that the first and third term on the right hand side 
of this expression cancel each other because of the identities 
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the gravitational acceleration in radial direction becomes 
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The corresponding expressions for the gravitational acceler- 
ation in the two angular directions are easy to obtain since the 
spherical harmonics Y ,m are the only angular-dependent terms in 
Eq. (fT4l i. Therefore, we only need to consider the partial deriva- 
tives of the spherical harmonics with respect to the 6 and co- 
ordinates. As the spherical harmonics are given by 



Y lm (0, cf>) = N lm P lm (cos 6)e im,p , 



(22) 



where N ,m is the normalization constant and P lm the associated 
Legendre polynomial, one finds 
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The derivatives of the associated Legendre polynomials are eas- 
ily obtained using the recurrence formula 



(x 2 - l)—P"\x) = lxP'"(x) - (I + m)P'"(x) . (25) 
dx 

Thus, one finds for the gravitational acceleration in the angular 
directions 
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^C lm (r) + r'D h "(r)\ (26) 



Annop Wongwathanarat et al.: An axis-free overset grid in spherical polar coordinates for simulating 3D self-gravitating flows 



5 



and 



1 



r sin 9 dip 



<D(r, 6, 



Y^-Y i,nY»» 



rsinO 



(0- 



—C' m (r) + r*D m (r)\. (27) 



Obviously, the expressions for three components of the gravita- 
tional acceleration (see Eqs. (f2TT >. d26b , and d27l >) are similar to 
that for the gravitational potential itself (see Eq. (TBI). Hence, 
besides computing derivatives of Legendre polynomials, our ex- 
tended Poisson solver can provide without much additional ef- 
fort both the gravitational potential and the corresponding accel- 
eration. 

Usage of the analytic expressions for the gravitational ac- 
celeration avoids the errors arising from the numerical differen- 
tiation of the gravitational potential. However, tests show that 
the results obtained using either the gravitational potential com- 
puted with the "standard" Poisson solver and subsequent numer- 
ical differentiation or directly the gravitational acceleration pro- 
vided by the extended Poisson solver differ only very slightly 
(see next section). Thus, we decided to stick to the "standard" 
Poisson solver in our simulations and compute the gravitational 
acceleration by numerical differentiation, as it requires no mod- 
ification of our code. 



4. Test Suites 

4.1. Sod Shock Tube 

The first problem of our test suite is the planar Sod shock tube 
problem, a classical hydrodynamic test problem (Sod 1978). We 
simulated this (ID Cartesian) flow problem using spherical co- 
ordinates and the Yin- Yang grid. The initial state consists of two 
constant states given by 



(p, P, v x ) 



(1.0,1.0,0.0) if x {n) > 0.4 cm 
(0.125,0.1,0.0) if x {n) < 0.4 cm 



(28) 



where p, p, and v x are the density, pressure and the velocity in 
x-direction of the fluid, respectively. We assume the fluid to obey 
an ideal gas equation of state with an adiabatic index y = 1 .4. 
The surface separating the two constant states is a plane orthog- 
onal to the x-axis located at jcW = 0.4 cm, the (positive) x-axis 
corresponding to a radial ray with angular coordinates 9^ = n/2 
and = 0. Thus, this ID planar Sod shock tube problem in- 
vokes all three spherical velocity components v r , vg, and when 
simulating the flow in spherical polar coordinates. This allows 
us to test both the scalar and vector transformations as well as 
the interpolation between the Yin and Yang grid patches. The 
simulation was carried out on an equidistant Yin- Yang grid of 
400 (r) x 92 (6) x 272 (0) x 2 zones (i.e., with an angular res- 
olution of one degree; see Eq.©). In radial direction the com- 
putational domain ranges from r = 0.05 cm to r — 1.0 cm. We 
impose a zero-gradient boundary condition at both edges of the 
radial domain. 

The solution of the shock tube problem is well-known. We 
compare our results with the so lution calculated using an ex- 
act Riemann solver (Toro 119971) . For comparison, data are re- 
sampled along the x-direction with a spacing Ax = 0.002 cm. 
Fig. [3] shows one dimensional profiles of p, p, v x , and e (spe- 
cific internal energy), respectively, along the x-direction at z (n) = 
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Fig. 3. One dimensional profiles of density p, pressure p, veloc- 
ity in x-direction v x , and specific internal energy e are shown 
along the x-direction at = 0.25 cm and = cm (dashed- 
dotted line in Fig.|4j) for the shock tube simulation at every 0.1 s. 
Open and filled symbols represent data points on the Yin and 
Yang grid, respectively. Solid lines give the distributions calcu- 
lated with an exact Riemann solver. 



0.25 cm and y (n) = cm (dashed-dotted line in Fig.[4]i at dif- 
ferent times. Our results agree very well with the solution ob- 
tained with the exact Riemann solver. The grid resolution is suf- 
ficiently high to give a sharp shock front and contact discon- 
tinuity while the rarefaction wave is smooth. The shock posi- 
tion is correct at all time throughout the simulation. The re- 
sampled data yield an accuracy of approximately 6% on average 
for shock positions. The shock wave and the contact disconti- 
nuity propagate smoothly across the Yin- Yang boundary located 
at x (n) = 0.25 cm without any noticeable effect by the existence 
of the boundary. To illustrate this behavior, Fig. [4] shows lines 
of constant density in the meridional plane (n) = at time 
t = 0.15 s. The isocontours are nearly perfectly straight lines 
perpendicular to the x-axis that are unaffected by the Yin- Yang 
boundary (dashed line). The contour lines are slightly bent near 
the outer radial edge of the computational domain due to the 
zero-gradient boundary condition we have imposed there. 

In order to firmly demonstrate that the Yin- Yang bound- 
ary does not cause numerical artifacts, we also computed this 
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Fig. 4. Snapshot of density contours in the meridional plane 
0W = at t = 0.15 s for the shock tube test problem. Dashed 
lines mark the Yin- Yang grid boundary, while the dotted circu- 
lar curves represent the inner and outer radial boundary of the 
computational domain, respectively. The one dimensional pro- 
files shown in Fig.[3]are re-sampled along the dashed-dotted line 
at z (n) = 0.25 cm. 



shock tube problem with a standard spherical polar grid using 
the same radial and angular resolution as for the Yin- Yang grid 
described above, i.e., 400 (r) x 180(0) x 360 (<p). We imposed re- 
flecting boundary conditions in 0-direction and periodic ones in 
0-direction. Fig. [5] shows a comparison of the results obtained 
with both simulations. The two panels give the tangential veloc- 
ity, defined as ^(y^) 2 + (vi"') 2 , in the meridional plane (n) = 
at time t = 0.15 s for the Yin- Yang grid (left), and the standard 
spherical polar grid (right), respectively. This velocity compo- 
nent should remain exactly zero because of the chosen initial 
conditions. Thus, it is a sensitive indicator whether the Yin- Yang 
boundary works properly, which obviously is indeed the case 
as the left panel of Fig. shows no hint of the location of that 
boundary. The modulus of the tangential velocity does nowhere 
exceed a value of 0.05 cm/s or approximately 5% of the shock 
velocity (in x-direction) except near the outer radial edge of the 
grids, where the boundary condition causes larger numerical er- 
rors. Note that nonzero tangential velocities are encountered on 
both the Yin- Yang grid and the standard spherical polar grid in 
the same grid regions at the same level. We thus conclude that 
they are the result of numerical errors that unavoidably occur 
when propagating a planar shock across a spherical polar grid, 
be it a standard one or a Yin- Yang grid. 

4.2. Taylor-Sedov Explosion 

As a second test for our code we consider the Taylor-Sedov 
explosion problem. We set up the initial state for the prob- 
lem by mapping a sph erically symmetric analytic solution 
dLandau & Lifshitz] 1 19591) onto the computational grid. We 
choose the parameters of the problem to mimic a supernova ex- 
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Fig. 6. Distributions of density (top), pressure (middle) and ra- 
dial velocity (bottom) versus radius from the explosion cen- 
ter (located at (x ( -"\y ( -"\z (n) ) = (7.0,0.0,2.5) x 10 19 cm for the 
Taylor-Sedov explosion problem plotted at every 10 11 s. Open 
symbols are data points from the Yin grid, while filled symbols 
represent sampled data from the Yang grid. The solid lines give 
the corresponding analytic solution. The data are re-sampled 
along the dashed-dotted line shown in Fig. [7] 

plosion in an interstellar medium. Because the shock wave re- 
sulting from the explosion is spherically symmetric with respect 
to the center of the explosion, we assume the explosion center to 
be located at the point {x {n \y {n) ,z (n) ) = (7.0,0.0,2.5) x 10 19 cm. 
Hence, this second test problem also involves a non-zero flux 
of mass, momentum, and energy across the Yin- Yang bound- 
ary, and as the previous shock tube test, it probes whether that 
boundary causes any numerical artifacts. 

The initial shock radius is ro = 2.9625 x 10 19 cm orrespond- 
ing to a time t exp = 0.34 x 10 11 s past the onset of the explosion, 
and the explosion energy was set to Eo = 10 51 erg. The ambient 
medium into which the shock wave is propagating is at rest. It 
has a constant density p/, = 10~ 25 g/cm 3 , and a constant pressure 
Pb — 1.4 x 10~ 13 erg/cm 3 . The fluid is described by an ideal gas 
equation of state with an adiabatic index y = 5/3, resulting in a 
density jump across the shock front of (y + l)/(y - 1) = 4. We 
use a grid resolution of 400 x 92 x 272 x 2 zones, a computational 
domain covering the radial interval r = [0.5, 15.] x 10 19 cm, and 
employ a zero-gradient boundary condition at both the inner and 
the outer radial boundary. 

Our results are shown together with the analytic solution in 
Fig. [6] We have re-sampled our data and calculated radial pro- 
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Fig. 5. Color maps of the tangential velocity defined by yj(v^) 2 + (v^) 2 in the meridional plane = resulting from the (ID 
Cartesian) shock tube problem. The snapshots are computed using the Yin- Yang grid (left) and a standard spherical polar grid (right) 
at a time t = 0.15 s. On the left panel, red and blue lines mark the boundaries of the Yin and the Yang grid patches, respectively. 
On the right panel, the two red circles show the inner and outer boundary in the radial direction of the standard spherical polar grid. 
The labels at the color bars give the tangential velocity in units of cm/s. The color range is limited to 0.05 cm/s to emphasize the 
smallness of the tangential velocity far from the outer radial grid boundary. 
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Fig. 7. Lines of constant density in the meridional plane </> (n) = 
obtained from our simulation of a Taylor-Sedov explosion. The 
snapshot is taken at a simulation time f J)m = 2.0 x 10 11 s which 
corresponds to an explosion time t exp ~ 2.34x 10 1 1 s. The dashed 
lines mark the Yin- Yang boundary, while the two dotted circles 
represent the inner and outer radial boundary of the computa- 
tional domain, respectively. The data presented in Fig. [6] are re- 
sampled along the dashed-dotted line . 



files of the density p, pressure p, and radial velocity v r along 
a line in z-direction through the explosion center using a uni- 
form radial spacing Ar = 10 18 cm. As one can see the numerical 
results agree very well with the analytic solution. All flow quan- 
tities are smooth across the Yin- Yang boundary, i.e., the shock 
wave passes that boundary without any noticeable numerical ar- 
tifact. Due to the finite resolution the density jump across the 
shock front is slightly smaller in the simulation than the analytic 
value of four. However, the shock front is sharp throughout the 
whole simulation, and it propagates with the correct speed. One 
distinct feature of the Taylor-Sedov solution is its spherical sym- 
metry. To illustrate that the Yin- Yang grid does not destroy this 
symmetry of the solution, we show a set of lines of constant den- 
sity in the meridional plane <n) = in Fig. [7] We also marked 
the line (dashed-dotted) along which the data given in Fig. [6] 
are re-sampled. The contour lines, all of which are almost per- 
fectly circular, are drawn at a simulation time f Jlm = 2.0 x 10 11 s 
(i.e., time step number 1276) corresponding to an explosion time 
f„ p ~2.34xl0 11 s. 

We further studied how the solution differs in the region 
where the Yin and Yang grid overlap. To this end we compare 
the total mass within the overlap region computed on the Yin 
and the Yang grid, respectively. Fig. [8] shows the evolution of 
the relative mass difference, i.e., the mass within the overlap 
region computed on the Yin grid, M°^ lp minus the mass com- 
puted on the Yang grid, M° vlp , divided by the sum of these two 
masses. We calculated this quantity for three different (angu- 
lar) grids with 400 x 32 x 92 x 2 zones (i.e., 3° angular reso- 
lution), 400 x 92 x 272 x 2 zones (i.e., 1° angular resolution), and 
400x182x542x2 (i.e., 0.5° angular resolution), respectively. For 
all three grid resolutions the relative mass difference has a value 
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Fig. 8. Evolution of the mass within the overlap region for the 
Taylor-Sedov test case computed on the Yin grid, M°^ lp minus 
the mass computed on the Yang grid, M° vlp , divided by the sum 
of these two masses. The dashed, dotted and solid lines give the 
solutions computed on a grid of 400 x 32 x 92 x 2 zones, (i.e., 3° 
angular resolution), 400 x 92 x 272 X 2 zones (i.e., 1° angular 
resolution), and 400 x 182 x 542 x 2 zones (i.e., 0.5° angular 
resolution), respectively. 



of about 10~ 4 . Although its evolution with time is different in 
case of the 3° simulation (because the coarse angular grid causes 
large errors when mapping the analytic initial data onto the grid 
which determine the further evolution), Fig. [8] shows that for an 
angular resolution better than 1° the relative mass difference be- 
haves similarly, its maximum value decreasing from 2.1 x 10~ 4 
at 1° angular resolution to 1.5 X 10~ 4 at 0.5° angular resolution. 

4.3. Rayleigh-Taylor Instability 

We also simulated a single mode Rayleigh-Taylor instability 
(RTI) on a Yin- Yang sphere. The initial configuration consists of 
a spherical shell of a heavier fluid of density p# = 2 g/cm 3 that 
is supported against a constant gravitational field g = 1 cm/s 2 
pointing in negative radial direction by a spherical shell of a 
lighter fluid of density pi — 1 g/cm 3 . The boundary between 
the two fluid shells is initially located at a radius r = 0.5 cm. To 
balance the gravitational force, the initial (radial) pressure dis- 
tribution is set to 



P(r) = 



P + gp H (l.0-r) 
P(r = 0.5) + gp L (0.5 ■ 



if r > 0.5 cm 
r) if r < 0.5 cm 



(29) 



where Pq — 1 erg/cm 3 . A radial velocity varying in angular di- 
rection as the spherical harmonics Y™(9, </>) with I = 3 and m - 2 
is used to perturb the initial configuration. The amplitude of 
the velocity perturbation is 2.5% of the local sound speed c s (r). 
Hence, the initial radial velocity is given by 



v r (r, 



-0.025 x c s (r) Yl(8, 0). 



(30) 



The spherical harmonics Y"'(9, </>) are connected with the associ- 
ated Legendre polynomials P"' via the expression 



Y?(6, 



m)\ 



(/ 



P'"(cos9)e im ^ . 



(31) 



t (s) 



Fig. 10. Position of the heads of the RTI bubbles versus time. 
Red symbols (circles, triangles, and squares) show data from the 
Yin grid, while blue symbols (diamonds) represent data on the 
Yang grid. 



The perturbation mode (/, m) 
velocity in the directions 



(3, 2) yields a maximum radial 



(9, <p) = {(tt - a, 0), (n - a, n), (a, n/2), (a, -n/2)} , (32) 

where a = arccos( V3/3). The remaining two velocity compo- 
nents of the perturbation mode are set equal to 0. The fluids are 
described by an ideal gas equation of state with an adiabatic in- 
dex y = 1.4. The simulation is carried out on a Yin- Yang grid 
of 400 x 92 x 272 x 2 zones. To keep the fluid in hydrostatic 
equilibrium, a zero-gradient boundary condition is used for both 
the inner and outer boundary in radial direction. The inner radial 
boundary is located at r — 0.1 cm. 

A snapshot of the resulting density distribution obtained with 
the Yin- Yang grid is displayed in Fig. [9] at epoch t = 2.85 s. The 
left panel shows color coded contour lines in 3D, and the right 
one a meridional cut at </> ( " ) = 0. The contour lines are drawn us- 
ing different color tables for the Yin and Yang grid, respectively. 
Four distinct bubbles of rising low density fluid (Yin: blue; Yang: 
bright gray) are clearly visible that reflect the initial perturbation 
mode (I, m) = (3, 2). High density fluid (Yin: yellow/red; Yang: 
dark gray/black) sinks down and settles at the inner part of the 
sphere. One can also notice Kelvin-Helmholtz instabilities de- 
veloping at the surface of the bubbles. This is particularly obvi- 
ous in the meridional cut (right panel). One of the RTI bubbles 
is within the Yang grid, while the three others reside on the Yin 
grid. It is obvious that the bubbles are distributed symmetrically 
following the perturbation pattern regardless of the grid patch. 
The 2D contour lines shown in the right panel of Fig.|9]empha- 
size this fact. 

The RTI bubbles grow with nearly the same growth rate in 
all four (perturbation) directions, as can also be seen from Fig. 
[10]that displays the position of each bubble's head versus time. 
The four curves lie exactly on top of each other during the phase 
of linear growth. There are slight discrepancies between the four 
curves in the non-linear regime, because the linear grid resolu- 
tion in angular direction is slightly non-equidistant (due to its 
9 dependence). Two curves from the Yin grid coincide perfectly 
since they represent the two bubbles that lie symmetrically above 
and below the equator in the Yin grid. The results confirm that 
the Yin- Yang grid does not favor any angular direction on the 
sphere. Since our aim was only to demonstrate this important 
fact, we do not further analyze the growth rate of the RTI. 
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Fig. 9. Surfaces of constant density in 3D (left) and 2D (right; meridional cut at = 0) resulting from the simulation of the 
Rayleigh-Taylor instability described in the text at t — 2.85 s. Contour lines on the Yin grid are shown using the blue-yellow colors 
while contour lines on the Yang grid are displayed using the white-black colors. 
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Table 1. Mean errors in the gravitational acceleration. 



4.4. Gravitational Potential of Homogeneous Spheroids 

We investigate the accuracy of our gravity solver by calculating 
the gravitational potential of homogeneous spheroids. We con- 
sider two homogeneous self-gravitating configurations: a pro- 
late spheroid with an axis ratio of 0.7, and a sphere. The con- 
figurations have a constant density p — 1 g/cm 3 , and are em- 
bedded into a homogeneous background of much lower density 
Pb = 10~ 20 g/cm 3 in order to minimize the background's con- 
tribution to the gravitational potential. The semi-major axis of 
the spheroid aligns with the x-axis, while its center is placed at 
the origin of the Yin- Yang grid. To provide a more difficult test 
for our multipole based gravity solver, we shift the center of the 
sphere off the origin of the computational grid by more than one 
sphere radius. 

The analytical form of the gravitational potential for both 
type of configurations are known. The solution for the prolat e 
spheroid can be found in chapter 3 of IChan drasekhar (1969), 
and the sphere's potential can be easily calculated. Fig.fTTIshows 
contour lines of the gravitational potential for both cases in the 
meridional plane <p {n) = 0. The potential is calculated on a grid 
of 400 x 92 x 272 x 2 zones with L = 15, where L is the number 
of spherical harmonics taken into account (see section [3~!2l . The 
contour lines are smooth across the Yin- Yang boundary for both 
the prolate spheroid and the sphere. 

Concerning the convergence behavior of the solver, we con- 
sider various grid resolutions and a number of spherical har- 



monics ranging up to L = 25 for this convergence test. The 
grid resolutions used in the test are 400 x 92 x 272 x 2 zones, 
400 x 47 x 137 x 2 zones, 200 x 92 x 272 x 2 zones, and 
200 x 47 x 137 X 2 zones, respectively. The maximum and mean 
error of the gravitational potential are given as a function of L 
for both considered configurations in the middle and right pan- 
els of Fig. QT| respectively. Both errors show a convergence be- 
havior with higher grid resolution, and tend to saturate at large 
values of L. This behavio r is similar to what is described in 
Miiller & Stei nmetzl (119951) . In addition, for lower grid resolu- 
tion the accuracy saturates at a lower number of spherical har- 
monics compared to calculations with a higher grid resolution. 
This is expected since higher order terms in the multipole expan- 
sion are not well represented on grids of lower angular resolu- 
tion. 

We also tested our extended Poisson solver discussed in sec- 
tion !3.2l In TableQ] we compare the mean errors in the com- 
ponents of the gravitational acceleration for both the prolate 
spheroid and the sphere test case computed with the numerically 
differentiated gravitational potential given in Eq. (fT4l with those 
obtained from the analytic expression given in Eqs. (ETT i. (|26| |. 
and d27j, respectively. We used a grid of 400 x 92 x 272 x 2 
zones and L — 15 for this comparison. 

For the prolate spheroid test case the "analytically" obtained 
accelerations exhibit a smaller mean error, especially for the 8- 
and ^-component of the gravitational acceleration. This results 
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Fig. 11. Contour lines of the gravitational potential (left column) for two homogeneous self-gravitating configurations: a prolate 
spheroid (top row) with an axis ratio of 0.7, and a sphere (bottom row). The configurations are indicated by the dark-gray shaded 
areas. Dashed lines show the Yin- Yang boundary, while dotted lines indicate the outer radial boundary of the computational grid. 
The middle and right columns give the maximum and mean error of the numerically calculated gravitational potential for different 
grid resolutions as a function of the number of spherical harmonics used in our multipole gravity solver. The solid, dotted, dashed, 
and dashed-dotted lines in both columns correspond to a grid resolution of 400 x 92 x 272 x 2 zones , 400 X 47 x 137 x 2 zones, 
200 x 92 x 272 x 2 zones, and 200 x 47 x 137 x 2 zones, respectively. 



from a strong decrease of the maximum error, which is large in 
regions where the angular components of the gravitational ac- 
celeration approach zero, i.e., near the major and minor axes of 
the prolate spheroid. However, in these regions the accelerations 
in 6 and ^-direction are orders of magnitude smaller than the 
radial component. Thus, they contribute only a tiny fraction to 
the total acceleration. In the sphere test case both variants of 
the extended Poisson solver produce similar mean errors. Based 
on these results we conclude that the extended Poisson solver, 
which provides the gravitational acceleration using analytic ex- 
pressions, works properly. Moreover, it gives a slightly more ac- 
curate gravitational acceleration, as it does not involve numer- 
ically differencing the gravitational potential. Nevertheless, for 
the reas ons stated in section[3~!2l we prefer to use the Poisson 
solver o f Mu ller&Steinmetzl (ll995 ) in our simulations. 

4.5. Self-gravitating Polytropes 

Using our Yin- Yang grid based hydro-code we have also con- 
sidered self-gravitating, non-rotating and rotating equilibrium 
polytropes. Both kinds of polytropes provide another test of the 
Poisson solver, and a test of how well our hydrodynamics code 
can keep a self-gravitating configuration in hydrostatic and sta- 



tionary equilibrium, respectively. In addition, the rotating poly- 
trope also serves to test the proper working of the Yin- Yang 
boundary treatment, as it involves a considerable and systematic 
flow of mass, momentum and energy flux across that boundary 
due to the polytrope's rotation. 

The polytropes have a polytropic index n — 1, a polytropic 
constant k = 1.455 x 10 5 , and a central density of p c = 7.905 x 
10 14 g/cm 3 . For our test runs we inte rpolated equ ilibrium poly- 
tropes calculated with the method of Eriguchi & Miiller ( 1985) 
onto a Yin- Yang sphere, and simulated their dynamic evolution 
(occurring as the interpolated configuration is not in perfect hy- 
drostatic equilibrium). The central region (r < 1 km) of the poly- 
trope is cut out and replaced by a corresponding point mass to 
allow for a larger time step. 

We use an artificial atmosphere technique to handle those re- 
gions of the computational grid that lie outside the (rotating, i.e., 
non-spherical) polytrope. The density in the atmosphere is set 
equal to a value p atm = 10~ 10 p o where p c is the central density 
of the polytrope. Here, atmosphere denotes any grid zone whose 
density is less than the cut-off density p cu t-off = 1C Pmax- 
Furthermore, for all zones in the atmosphere the velocity is set to 
zero in order to keep the atmosphere quiet. This procedure is ap- 
plied at the end of every time step throughout the simulation. A 
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Fig. 12. Relative change of the central density of a non-rotating 
(nearly) equilibrium polytrope as a function of time. 



zero-gradient boundary condition is imposed at the outer radial 
boundary, and a reflecting boundary condition at the inner one. 
The polytrope's evolution is followed for 10 ms corresponding 
to approximately 10 dynamic time scales in order to check how 
well the initial approximate equilibrium configuration is main- 
tained by the Yin- Yang code. 

For the non-rotating polytrope, we employ a grid of 400 x 
20 x 56 x 2 zones. Note that we are able to use a relatively low 
angular resolution compared to the other tests, because the prob- 
lem has spherical symmetry. Our results show that the polytrope 
stays perfectly spherically symmetric throughout the simulation, 
and that the non-radial velocities inside the polytrope remain 
zero. This demonstrates that the Yin- Yang grid is able to pre- 
serve the initial spherical symmetry. Fig.[T2]shows the evolution 
of the central density (more precisely of the density of the inner- 
most radial zone at r — 1 km), which exhibits oscillations with 
an amplitude of the order of 10~ 4 without any sign of a system- 
atic trend. Comparing the initial radial distributions of the den- 
sity (Fig. [13] upper panel) and the radial velocity (Fig. [13] lower 
panel) of the polytrope with those after 10 ms of evolution, we 
find no significant deviations. Relative changes in the density 
profile are of the order of 10~ 4 , comparable to the size of the 
fluctuations of the central density. Only for zones near the edge 
of the polytrope the deviations can reach a level of up to 20%, 
in particular in the zone next to the atmosphere. The figure also 
shows that data points from the Yin and the Yang grid lie on top 
of each other confirming that the code preserves the initial spher- 
ical symmetry of the polytrope very well. Except for the zones 
at the polytrope's surface, where the radial velocity is fluctuating 
at a level of approximately 2 x 10 8 cm/s, the radial velocities are 
less than 10 6 cm/s (i.e., less than 0.1% of the local sound speed). 
Thus, we conclude that a non-rotating (n = 1) equilibrium poly- 
trope is correctly handled by our Yin- Yang hydro-code. 

The rotating polytrope needs a higher grid resolution in 9- 
direction, as it is no longer spherically symmetric. Thus, we used 
a grid resolution of 400 x 92 x 272 x 2 zones for this simulation. 
The initial oblate equilibrium configuration has an axis ratio of 
0.7. We, again, evolve the configuration for 10 ms to test the cor- 
rect treatment of the situation by our Yin- Yang hydro-code. 

Fig. [14] shows the relative variation of the central density as 
a function of time along an equatorial ray (# (n) = tt/2; = 0) 
and along the pole = 0), respectively. One also recognizes a 
slight systematic trend in the behavior of the density fluctuation, 
which is steeper along the equator than at the pole. However, 
in both cases the relative increase of the central density is very 
small (~ 10~ 3 ). The initial radial density profiles along the pole 
and the equator do not show any significant change during the 
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Fig. 13. Density (top) and radial velocity (bottom) of a non- 
rotating n = 1 equilibrium polytrope as a function of radius after 
t = 10 ms of "evolution". In the top panel, the solid line shows 
the initial density profile. Red circles and blue triangles corre- 
spond to data from the Yin and the Yang grid, respectively. 
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Fig. 14. Same as Fig. [12] but for a rotating polytrope. The solid 
and dashed curves show the relative variation of the density 
along an equatorial ray (# (n) = 7r/2; </> ( "- ) = 0) and along the 
pole (6* (n) = 0), respectively. 



10 ms of evolution we have simulated with the Yin- Yang code 
(Fig. US upper panel). The axis ratio has slightly increased to 
a value of 0.719. The radial velocities (Fig. [15] lower panel) 
are larger than in the non-rotating case by about an order of 
magnitude, because it is obviously more difficult to keep a ro- 
tating polytrope in equilibrium than a non-rotating (spherically 
symmetric) one. We again find the largest radial velocities (a 
few times 10 8 cm/s) near the surface of the polytrope, especially 
along the equator. However, these velocities vary with time. 
When averaged over time (in the time interval t = [9, 10] ms) the 
profiles become flatter and the velocities smaller. This confirms 
that the polytrope is oscillating around its equilibrium configu- 
ration. 
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Fig. 15. Density (upper panel) and radial velocity (lower panel) 
of a « = 1 rotating polytrope in stationary equilibrium as a func- 
tion of radius after t — 10 ms of "evolution". In both panels the 
solid and dashed lines show the profiles along an equatorial ray 
(0 (n) = tt/2, 4> {n) = 0) and along the pole (# (w) = 0), respec- 
tively. Red circles and blue triangles in the upper panel corre- 
spond to data from the Yin and the Yang grid, respectively. In 
the lower panel, we show in addition time averaged (over the in- 
terval t = [9, 10] ms) velocity profiles along the equatorial ray 
(dotted) and the pole (dashed-dotted). 



5. Conservation problem 

The Yin- Yang grid has a disadvan tage common with other types 
of overlapping grids (see, e.g., Che s shire & Hen shawl Il994t 
IWangll 19951; IWu et alJl2007D . The communication via interpola- 
tion between the two grid patches does not guarantee conserva- 
tion of conserved quantities even though the finite-volume differ- 
ence scheme employed on each grid patch is conservative. Non- 
conservation occurs when a flow across the Yin- Yang boundary 
is present. This is the case in most of our tests except for the 
simulation of the non-rotating polytrope that involves only ra- 
dial flow. 

Nevertheless, we are still able to obtain sufficiently good re- 
sults for all the test simulations discussed in the previous section. 
The degree of non-conservation is highly problem dependent. A 
simulation involving a considerable and systematic flow across 
the Yin- Yang boundary, as e.g., in the case of the rotating poly- 
trope, will result in a larger degree of non-conservation. We ob- 
serve that the total mass increases by 0.07% within 10 ms (or 
about ten dynamical timescales) in the case of the rotating poly- 
trope. For the Taylor-Sedov test case, which is the cleanest test 
case in this respect (as it involves, e.g., no boundary effects like 
the shock tube, and e.g., no artificial atmosphere like the rotating 
polytrope), we find a mass loss of the order of 10~ 5 , only. As 
Fig. Q2] demonstrates this mass loss can be reduced by using a 
higher angular resolution. 

Conservation of conserved scalar quantities can be obtained 
to machine precision by applying the algorithm described in de- 
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Fig. 16. Evolution of the relative mass loss, (M—Mq)/Mq, where 
Mo is the initial total mass, for the Taylor-Sedov test simulated 
on three different (angular) grids with 400 x 32 x 92 x 2 zones 
(i.e., 3° angular resolution; dashed line), 400x92x272x2 zones 
(i.e., 1° angular resolution; dotted line), and 400 x 182 x 542 x 2 
zones (i.e., 0.5° angular resolution; solid line), respectively. 



tail in iPeng et al.1 (12006). and summarized below. According to 
this algorithm scalar fluxes at the outer edges of boundary zones 
of both the Yin and the Yang grid are replaced by scalar fluxes 
computed using only "interior" fluxes from adjacent grid zones. 

As an illustration, consider the Yin- Yang grid overlap con- 
figuration in Fig. [17] where PQRS is a grid zone at the boundary 
of the Yang grid (blue) which overlaps with the underlying grid 
zone ABCD of the Yin grid (red). Fluxes referring to the Yin and 
the Yang grid are denoted by / and g, respectively. 

The boundary flux gpQ of the Yang grid is replaced by the 

flux 

fpQ = Sfq + fpF, (33) 

where fp q and fpp are the fluxes through the segments FQ and 
PF, respectively. 

The flux fp q in Eq.(l33l is calculated using information from 
zone ABCD. The evolution of a scalar quantity ^abcd of zone 
ABCD is given by 



£aBCD - %ABCD + (fAB ~ fcD + fBC ~ /ad)- 



(34) 



Similarly, for the fraction of the zone ABCD defined by the poly- 
gon ABFED one has, 



DE i>i 

%'abfed - Zabfed + (/ab~ fcD== +/bc= ~/ad -/ef)- (35) 



BF 
BC 



Assuming a piecewise constant state within the zone ABCD, 
Eqs.OD) and ([35]) lead to 

a ^ABCD ~ %ABCD> ~ £ ' ABFED ~ £ ' ABFED (^6) 

where a is the overlapping volume fraction (area) described in 
section 3. Therefore, 

DE W 

o(fAB - fcD + fBC - /ad) - fAB - fcD + /fiC^= _ /AD - fEF 

CD BC 

(37) 

Note that the flux fpp is the only unknown in Eq.OTb. Since the 
intersection points E and F are already known from the step to 
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Fig. 17. Illustration of the Yin- Yang grid overlap configuration, 
where PQRS is a grid zone at the boundary of the Yang grid 
(blue) which overlaps with the underlying grid zone ABCD of 
the Yin grid (red). Fluxes referring to the Yin and the Yang grid 
are denoted by / and g, respectively. 

calculate the volume fraction a, the lengths of all segments can 
be obtained. The flux fp q is then given by 



FQ 

f FQ -f EF =. 



(38) 



After obtaining the still missing flux fpp in Eq.(l33l by a similar 
procedure, the scalar quantity %pqrs of the boundary zone PQRS 
is updated according to 

= ?PQRS + (8QR ~ 8PS + gRS ~ f PQ )- (39) 

This procedure is then repeated to update all boundary grid 
zones. 

After implementing the above algorithm we are able to con- 
serve mass and total energy up to machine precision. However, 
the conservation of momentum is more complicated since the 
momentum equations in spherical coordinates involve not only 
flux (i.e., divergence) terms but also source terms (due to the 
presence of fictitious and pressure forces), and due to the "mix- 
ing" of momentum components as the Yin and Yang grid patches 
are rotated relative to each other (see Fig.[T7l>. 

As we have not yet devised and implemented a correspond- 
ing momentum conservation algorithm, momentum is not yet 
perfectly conserved in our code. For that reason we also refrain 
from using the scalar conservation algorithm described above, 
since in some simulations {e.g., in the Taylor-Sedov explosion 
simulation) we encountered a negative internal energy in some 
zones due to the inconsistency arising from the perfect conser- 
vation of mass and total energy on one hand and the imperfect 
conservation of momentum on the other hand. In our test runs the 
momentum violation is small, e.g., amounting to 0.24% (0.03%) 
angular momentum loss in the case of the rotating polytrope for 
a grid with three (one) degree angular resolution. 

6. Performance and Efficiency 

One of the main purposes in implementing the Yin- Yang grid is 
to ease the severe restriction imposed on the size of the time step 
for any explicit hydrodynamics scheme by the CFL condition in 
the polar regions of 3D simulations using a grid in spherical po- 
lar coordinates. In most applications the size of the time step is 



computational domain 


angular grid resolution 


gain factor 


tiill A.TT *;nhprp 

114.11 TVL auiicic 


3° 


26 




2° 


40 




1° 


80 


sphere except for a cone 
of 5° half opening angle 
cut-out at both poles 


1° 


7 



Table 2. Expected gain factor when using the Yin- Yang grid. 



restricted most strongly by the size of the zones in ^-direction, 
which is smaller than the size in ^-direction by the factor sin 8 as- 
suming an equal angular resolution 5 = AO — A<p in both angular 
directions. 

For a spherical polar grid the factor sin 9 implies (assuming 
zone centered variables) a minimum zone size 

df = 6 sin(«5/2) 

(in radians) in ^-direction for the first zone next to the pole. 
Typically, sin(5/2) « 10~ 2 . On the other hand, applying the Yin- 
Yang grid yields 



dl Y = 5 sin(7r/4 • 



512) 



for the size of the smallest zone in ^-direction, which is typically 
about 0.7. Hence, for the Yin- Yang grid the smallest zone size in 
azimuthal direction is larger by the ratio 



,sph 



sin(7r/4 - 5/2) 
sin(5/2) 



(40) 



compared to the spherical polar grid. 

Table[2]gives the value of this ratio for grids of various angu- 
lar resolution, and various computational domains. These num- 
bers provide an estimate of the gain in computation time one can 
expect when using the Yin- Yang grid instead of the spherical 
polar grid. 

However, the gain factor calculated from the relative grid 
spacings does not determine the gain in the size of the time step, 
as the latter is given in a more complicated way by the CFL con- 
dition 



At, 



CFL 



< C 



V(l 



rA9 



r sin 6A(p 



Ar 2 +(rA0) 2 +(rsin0A0) 2 



,-i/2 



(41) 



where C, v r , vq, v^, and c s are the Courant factor, the flow veloc- 
ities in radial, colatitude and azimuthal direction, and the local 
sound speed, respectively. The CFL condition shows that the in- 
crease in the size of the CFL time step is somewhat smaller than 
implied by the gain factor resulting from the ratio of the sizes of 
the smallest zones of the Yin- Yang grid and the spherical polar 
grid. In addition, the increase of the time step is problem depen- 
dent. 

Besides the performance gain due to the increased size of the 
CFL time step, the Yin- Yang grid also requires less computa- 
tional zones to cover the full sphere, and thus less computational 
time. For an angular resolution 6 the spherical polar grid needs 

(tt/6) x (2n/6) 

zones to cover the full sphere, while the Yin- Yang grid requires 
only 

(tt/2S + 2) x (3n/26 + 2) x 2 
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zones. Hence, up to 25% fewer computational zones are re- 
quired. The gain depends only weakly on angular resolution and 
is problem independent. 

However, employing the Yin- Yang grid also requires some 
extra amount of computation compared to the spherical polar 
grid (see Sec. 3). In the following we only consider the extra 
costs of calculations during the actual simulation, but not the 
extra costs arising during the initialization, since these are neg- 
ligible. We emphasize again that there are two major extra sets 
of calculations necessary when applying the Yin- Yang grid. The 
first set concerns the interpolation of the ghost zone values that 
are needed for the communication between the Yin and Yang 
grid patches. The second set arises from the interpolation of the 
density onto the auxiliary spherical polar grid grid and the in- 
terpolation of the gravitational potential back from the auxiliary 
grid onto the Yin- Yang grid. Exploiting the algorithms described 
in Sec. 3, the computational cost for both parts is almost negli- 
gible compared to the total computing time. Interpolation of the 
ghost zone values requires only 2.3% of the total computing time 
per cycle in simulations with self-gravity, while the interpolation 
of density and gravitational potential performed within the grav- 
ity solver accounts for 1.5% of the computing time needed for 
the gravity solver. This corresponds to approximately 0.3% of 
the computing time per cycle. 

To obtain actual numbers for the gain, we performed several 
timing tests including simulations with and without self-gravity 
using four different grid resolutions. The tests were carried on an 
IBM Power6 using a single processor. According to these tests 
the computing time per cycle for the Yin- Yang grid averaged 
over five cycles is approximately 15% and 20% smaller than for 
the spherical polar grid for simulations without self-gravity and 
with 2° and 1 angular resolution, respectively. For simulations 
including self-gravity, the gain factor decreases by 3% approxi- 
mately. 

Concerning the gain from the less restrictive CFL condition, 
we consider the case of the rotating polytrope since the size of 
the time step does not vary much throughout the simulation. For 
an angular resolution of 1°, we find a gain of approximately a 
factor of 63 when using the same Courant number both for the 
Yin- Yang grid and the spherical polar grid. 

7. Conclusion 

A two-patch overset grid in spherical coordinates called the 
"Yin- Yang" grid is successfully implemented into our 3D 
Eulerian explicit hydrodynamics code, PROMETHEUS, includ- 
ing in particular the treatment of self-gravitating flows. The Yin- 
Yang grid eases the severe restriction of the time step size in the 
polar regions of the sphere, because each Yin- Yang grid patch 
contains only the low-latitude part of the usual spherical polar 
grid. From our experiences, the implementation steps are easy 
and straightforward for a hydrodynamics code using directional 
splitting and having 3D spherical polar coordinates already im- 
plemented due to the simplicity of the Yin- Yang transformation 
and its symmetry property. Basically it involves doubling the 
state variable arrays, calling the ID core hydrodynamics solver 
in angular directions for both the Yin and the Yang grid, and 
adding a subroutine that handles the Yin- Yang transformation 
and the interpolation of variables between both grids. 

We validated the code with several standard hydrodynamic 
tests. The test results show good agreement with analytic so- 
lutions if these are available. Furthermore, as demonstrated by 
three of our test problems - a planar shock crossing the Yin- 
Yang grid boundary, an off-center Taylor-Sedov explosion in- 



volving mass, momentum (all three components) and energy flux 
across the Yin- Yang grid boundary, and a polytrope whose rota- 
tion leads to considerable and systematic mass, momentum (only 
angular components) and energy flux across the Yin- Yang grid 
boundary - the Yin- Yang grid does not introduce any numer- 
ical artifact at the internal Yin- Yang boundary. The tests also 
confirm that the numerical solutions obtained with the Yin- Yang 
grid do not show any evidence of a preferred radial direction, 
as it eliminates numerical axes artifacts which seriously flaw the 
flow near the coordinate symmetry axis when using a spherical 
polar grid. Besides successfully simulating a Taylor-Sedov ex- 
plosion and self-gravitating (rotating and non-rotating) equilib- 
rium polytropes the code has also passed another astrophysically 
relevant test involving the growth of Rayleigh-Taylor instabili- 
ties. 

Because the communication between the two grid patches in- 
volves interpolation, flows across the Yin- Yang boundary cause 
some small amount of non-conservation of conserved quanti- 
ties. However, even for the (in this respect) severe test case of 
the rotating polytrope involving large flows across the Yin- Yang 
boundary, we observe only a small amount (less than one per- 
cent) of non-conservation. 

The Yin- Yang grid offers a large gain in computing time 
arising from two sources. Firstly, the number of computational 
zones needed is reduced by 20% approximately depending on 
the angular resolution. This gain reduces the computing time per 
cycle and is problem independent. Secondly, the size of the CFL 
time step is considerably enhanced, because the polar regions 
with converging meridional coordinate lines are not present in 
case of the Yin- Yang grid. The corresponding gain in time step 
size highly depends on the problem simulated. The extra costs 
for interpolation between the two grid patches and the interpola- 
tion performed in the gravity solver are negligible compared to 
the gain in the time step size. 

In conclusion, our implementation of the Yin- Yang grid into 
the multi-dimensional hydrodynamics code PROMETHEUS 
brings about the possibility to simulate three dimensional self- 
gravitating hydrodynamic flows in spherical coordinates which, 
in most cases, have been computationally inaccessible up to now 
due to the prohibitively large computational costs. With the pos- 
sibility to add more physics such as neutrino transport (work 
in progress), the new code version can be used to carry out, 
e.g., core collapse supernova simulations in 3D. 
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